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Abstract 

A  numerical  technique  is  developed  for  solving  coupled  electrochemical  reaction-diffusion  equations.  Through  analyzing  the 
nonlinearity  of  the  problem,  a  trial  and  error  iterating  procedure  is  constructed.  The  coefficient  matrix  is  arranged  as  a  tridiagonal  form  with 
elements  of  block  matrix  and  is  decomposed  to  LU  form.  A  compact  forward  and  backward  substitution  algorithm  based  on  the  shift  of 
inversing  block  matrix  by  Gauss-Jordan  full  pivoting  is  developed.  A  large  number  of  node  points  is  required  to  converge  the  calculation. 
Computation  experiences  show  that  the  iteration  converges  very  quickly.  The  effects  of  inner  diffusion  on  the  electrochemical  reaction  are 
analyzed  by  numerical  solutions.  ©  2002  Published  by  Elsevier  Science  B.V. 
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1.  Introduction 

Industrial  chemical  reactions  are  usually  accompanied 
with  mass  and  energy  transfer,  either  homogeneously  or 
heterogeneously.  Mathematical  modeling  for  these  pro¬ 
cesses  is  based  on  material  and  energy  balance.  One  can 
generate  a  set  of  differential  equations  known  as  the  reac¬ 
tion-diffusion  problem.  Owing  to  the  strong  nonlinearity  of 
the  reaction  rate,  mainly  from  the  effect  of  temperature, 
reaction-diffusion  equations  are  paid  more  attention  in 
analyzing  and  designing  chemical  and  catalytic  reactors 
and  are  the  major  role  in  analyzing  the  nonlinear  dynamic 
behaviors  in  reactor  engineering.  The  same  phenomena  exist 
in  electrochemical  processes,  with  the  add  complexity  of  a 
varying  potential  field,  and  considerable  research  has  been 
reviewed  for  electrochemical  reactions  occurring  in  the 
porous  electrode  [1], 

Newman  discussed  the  numerical  solution  of  boundary- 
value  problems  consisting  coupled  ordinary  differential 
equations  which  one  can  often  met  in  chemical  engineering 
science,  and  developed  a  unique  technique  to  solve  coupled, 
linear  differential  equations  [2-4].  In  his  procedure,  a  large, 
sparse  matrix  was  collapsed  to  a  tridiagonal  matrix  with 
elements  of  block  matrixes,  which  make  it  easy  to  be 
inversed  in  solving  the  algebraic  equations  transformed  from 
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the  original  differential  equations.  Accordingly,  the  subrou¬ 
tine  in  his  numerical  method  is  called  Newman’s  Band(j). 
White  discussed  Newman’s  method,  and  promoted  this 
technique  for  wide  application  [5-7],  However,  for  solving 
reaction-diffusion  equations — the  nonlinear  two-point 
boundary  value  problem,  respective  trial  and  error  algo¬ 
rithms  are  needed.  Shooting  methods,  such  as  Runge-Kutta 
integration  scheme,  is  a  marching  technique  that  must  be 
backwards  processed  from  end  point  to  beginning  point  to 
avoid  inherent  instability.  This  method  is  commonly  used, 
but  it  depends  on  a  proper  choice  of  the  initial  value  at  the 
ending  point.  Linearizing  the  nonlinear  kinetics  term  that 
convert  the  nonlinear  differential  equations  into  a  linear  one 
is  called  the  linearization  trial-and-error  technique,  but  it 
does  not  guarantee  that  all  steady  state  solutions  can  be 
determined.  Compared  with  these  two  methods,  transform¬ 
ing  the  nonlinear  differential  equations  into  nonlinear  alge¬ 
braic  equations  that  are  then  solved  by  the  Newton  method 
did  not  receive  the  appropriate  attention  because  the  calcu¬ 
lating  time  and  storage  is  greater. 

In  this  paper,  we  extend  Newmans’s  matrix  method  to 
solve  coupled  electrochemical  reaction-diffusion  equa¬ 
tions — a  coupled  high  nonlinear  two-point  boundary  value 
problem  that  is  encountered  in  modeling  of  a  fuel  cell 
catalyst  layer.  Based  on  Taylor  expansion,  a  trial  and  error 
iteraction  algorithm  is  developed  for  solving  the  nonlinear 
algebraic  equations  transformed  from  the  electrochemical 
reaction-diffusion  equations.  Owing  to  the  adoption  of 
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Nomenclature 

a  specific  interfacial  area  (cm-1) 

Co  dimensionless  form  of  CH+  (co  =  CH+/C^J 

Ci  dimensionless  form  of  Ch2  (ci  =  Ch2/Ch2) 

Ct  dimensionless  electric  potential  in  the  solid 
matrix  phase  (C2  =  «/<"/> , ) 

C3  dimensionless  electric  potential  in  the  pore  ionic 
phase  (C3  =  nfcj) 2) 

CH+  effective  concentration  of  proton  per  unit  volume 
of  solution  (mol  cm-3) 

Ch2  effective  concentration  of  hydrogen  per  unit 
volume  of  solution  (mol  cm-3) 
initial  effective  concentration  of  hydrogen 
(mol  cm-3) 

D,l2  effective  diffusion  coefficient  of  hydrogen  in  pore 
phase  (cm2  s-1) 

E°  formal  potential  (V) 

/  constant  (/  =  F/RT  (V-1)) 

F  Faraday’s  constant  (96485  C  eq.-1) 

i  total  current  density  leaving  the  matrix  phase 

(A  cm-2) 

z'i  superficial  electronic  current  density  in  the 
matrix  (A  cm-2) 

i2  superficial  ionic  current  density  in  pore  phase 
(A  cuT2) 

/ceii  current  density  of  fuel  cell,  defined  as  positive 
(A  cm-2) 

k°  heterogeneous  rate  constant  for  hydrogen  oxidation 
(cm  s-1) 

k\  dimensionless  parameter  (k\  =  nflIce\\/o) 

k2  dimensionless  parameter  (k2  =  nflIceu/K) 

I  thickness  of  anode  catalyst  layer  (cm) 
n  =  2  number  of  electrons  transferred  in  electrode 
reaction  (=2) 

,V|[,  superficial  flux  density  of  hydrogen 
(mol  cm-2  s-1) 

R  universal  gas  constant  (8.314  J  mol-1  K-1) 

RHl  electrode  reaction  rate  of  hydrogen,  defined  as 
positive  (mol  cm-3  s-1) 

T  absolute  temperature  (K) 

U  dimensionless  formal  potential  ( U  =  nfE°  ) 
x  dimensionless  form  of  distance  through  porous 
electrode  (x  =  X/I) 

X  distance  through  porous  electrode  (cm) 

Greek  symbols 
a.  transfer  coefficient 

cf) i  electric  potential  in  the  solid  matrix  phase  (V) 

(/> 2  electric  potential  in  the  pore  ionic  phase  (V) 

<P  Thiele  modulus  (Z>2  =  cil2k° /DUl) 

<F\  dimensionless  parameter  (<P2  =  an2fFk0l2Cft2/(T) 

<k>2  dimensionless  parameter  =  an2fFk°l2C^2/K) 

k  effective  conductivity  of  the  pore  ionic  (proton) 

phase  (Q-1  cm-1) 

er  effective  conductivity  ofthe  solid  matrix  (O  '1  cm-1) 


Newman’s  unique  method,  in  each  iteration  loop  a  compact 
forward  substitution  can  be  processed.  In  each  forward  or 
backward  substitution,  we  use  Gauss-Jordan  full  pivoting 
technique  to  guarantee  the  numerical  stability.  The  whole 
computing  is  very  fast,  even  though  we  set  2000  node  points 
to  make  the  substitution  effectively.  For  a  general  case  of 
an  electrochemical  reaction-diffusion  process  in  a  fuel  cell 
catalyst  layer,  just  two  iteration  loops  can  converge  the 
solution.  Using  this  method  to  study  the  effect  of  Thiele 
modulus  on  the  electrochemical  reaction-diffusion  process, 
some  intrinsic  phenomena  due  to  the  nonlinear  behavior  of 
the  system  is  found  by  numerical  solution  that  is  beneficial  to 
analyze  and  scale-up  a  fuel  cell  system. 


2.  Model  equation 


We  consider  for  the  isothermal  hydrogen  oxidation  reac¬ 
tion  occurring  in  a  porous  catalyst  layer  of  a  fuel  cell  anode. 

H2^2H+  +  2e-  (1) 


A  schematic  representation  of  it  is  shown  in  Fig.  1.  The 
catalyst  layer  is  viewed  as  a  continuum  of  two  phases,  each 
phase  either  a  pure  ionic  or  electronic  conductor.  The  ionic 
phase  is  considered  a  cationic  selective  polymer,  such  as 
Nafion.  Therefore,  the  proton  concentration  throughout  the 
catalyst  layer  is  fixed.  This  uniform  electrolyte  concentra¬ 
tion  means  the  superficial  current  density  in  the  pore  ionic 
phase  is  due  only  to  ion  migration.  This  is  represented 
mathematically  as 


h  =  -k 


#2 

dX 


(2) 


The  rate  at  which  ionic  current  enters  the  pore  solution,  di2l 
dX,  is  proportional  to  the  reaction  rate  on  a  volumetric  basis 
that  is  expressed  by  Butler- Volmer  expression  (defining 
anodic  current  as  positive)  [8] 

=  nFRn2  =  anFk°\Cu1e('-a)nf(ll,'-rl,2-!'j/) 

-  CH+e-m/Wl-^-£0,)]  (3) 


The  movement  of  electrons  in  the  solid  matrix  phase  of  the 
porous  electrode  is  governed  by  Ohm’s  law 


i\  =  —a 


dX 


(4) 


The  consequence  of  electroneutrality  is  that  the  divergence 
of  the  total  current  density  is  zero. 


d;  dzj  d/'2 
dX~  dX  +  dX 


(5) 


The  flux  of  dissolved  hydrogen  in  the  porous  anode  is 
determined  by  diffusion 


^h2 


—  —Dh2 


dCn2 

dX 


(6) 
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Fig.  1 .  Hydrogen  oxidation  in  a  fuel  cell  anode. 


Material  balance  for  this  reaction-diffusion  process  can  be 
modeled  according  to  equation  of  continuity  at  steady  state  [9] 


cPVh, 

dX 


=  Ru2 


Substituting  Eq.  (3),  (6)  into  Eq.  (7)  yields 

dh2  =  ak°[CH2e(l-a)nf(^--By) 


(7) 


(8) 


Differentiating  Eqs.  (2)  and  (4)  and  substituting  Eq.  (3)  into 
the  resulting  equation  gives 

o 

dXz 

-  (9) 

Kpfy  =  -nFakQ[Cvl2e(l-(l)nf(^-,l,^) 

dX1 

-  (10) 


Boundary  conditions: 

1.  @  X  =  0, 


oS 

II 

(11) 

£  j1 

II 

1 

Q  In'1 

(12) 

<Pl  =  0 

(13) 

2.  @  X=l, 

ii 

O 

(14) 

o 

II 

(15) 

d(f>2  ^cell 

dX  k 

(16) 

Setting  the  following  dimensionless  variables  and  parameters 

Ch+  Ch2  ,, 

Co  =  ~^T ,  Cl  =  -Q- ,  c2  =  nf4>u  c3  =  nf<p2, 

lh2  lh2 


X  9  kal2 

x  =  —,  <P2  =  a - (Thiele  modulus), 

l  Du, 

b-0 12  r'O 


U  =  nfEr 


<k\  = 


an2fFk°!2C^ 


2  anzfFk'rC'u2 


n  = 


,  nfl 
*i  =  — rcen, 
er 


,  nfl 

ki  =  —  fee  n 

K 


Then  Eqs.  (8),  (10)  and  (11)  are  reduced  to  the  dimensionless 
form 

^  =  ^2[c1e(1-a)(e2-C3-t/)  -  coe~a{c2-C3~u)]  (17) 

axz 

=  <2>2[CieU^)C2-C3-E/)  _  Coe- a(c2-C3-t/)]  (18) 

^  =  -&[c1eQ-aX‘*-,*-W  -  c0e~a^-C2-u)]  (19) 

d\” 

Boundary  conditions: 


1 .  @  x  =  0, 

Cl  =  1 

(20) 

dc2  , 

d*  =  ~k' 

(21) 

C3  =  0 

(22) 

2.  @  x  =  1, 

o 

II 

•8  l-s 

(23) 

dC2 

df=° 

(24) 

1 

II 

SI-3 

(25) 

Although  the  above  model  Eqs.  ( 17)— (25)  are  for  an 
isothermal  system,  due  to  the  unique  nonlinear  character¬ 
istics  of  the  electrochemical  reaction  term  that  is  similar 
to  the  effect  of  temperature  on  a  chemical  reaction  rate, 
the  numerical  solving  method  for  these  electrochemical 
reaction-diffusion  equations  can  be  transformed  to  study 
non-isothermal  problems. 
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3.  Numerical  solution 

For  solving  the  model  Eqs.  (17)-(25),  they  are  cast  into 
sets  of  nonlinear  algebraic  equations  approximated  by 
finite  difference.  Then  an  iterative  algorithm  is  constructed 
based  on  the  Taylor  expansion.  However,  if  one  adds  the 
Eqs.  ( 1 7)— ( 1 9)  together  to  try  to  eliminate  the  nonlinear 
term,  the  computation  is  not  stable.  Analyzing  the  com¬ 
puting  process  on  a  microcomputer  with  eight  word  bits 
for  a  double  precision  real  showed  that  the  block 
matrixes  were  ill-conditioned.  The  condition  numbers 
are  too  large  to  inverse  the  coefficient  block  matrixes 
correctly  due  to  the  error  accumulation  even  for  the  second 
node  point. 

Therefore,  at  each  node  point  except  the  boundary  points, 
the  second-order  derivatives  can  be  approximated  by  three 
point  central  difference  accurate  to  0(h2).  For  example, 
Eq.  (17)  is  approximated  as 

ci  O'  -  1)  -  2d  O')  +  CiO+  !) 

-  A2^2(ci(y>(1_“)[c2W_C3W_t/1-  Coe-a[c2(J)-C3(i)-u])  =  0 

(26) 

If  there  are  enough  node  points,  we  can  suppose  that  the 
difference  equation  at  node  point  j,  such  as  Eq.  (26),  depends 


only  on  the  variables  at  three  points.  If  the  left-hand  term 

is  expressed  as  F\  .-,  then  the  Taylor  series  expands  it  to 

2 

accuracy  0(h  )  as 


=  c?(;'-l)-2c?(j)  +  c°(;+l) 

-  h2<p2[c°l(j)e(l^a)^('j)~c°^j',~u]  -  coc~a[c“(7)_e°(7)_£/|] 

+  [ct0-l)-c?0-l)]  dFl,i 
+  [d0+i)-c?0+i)] 

+  foO)  -  c°0)] 


dc  1,7-1 

dFij 


.Of  ^dFU 


dcXj 


dc\,j+\ 
=  0 


+hU)-c?U)& 

o  ochi 

+  [C2U)-C°2(j)}p^ 

0  CfC2  J 


(27) 


Substituting  corresponding  partial  derivatives  and  intro¬ 
ducing 


ro  =  e  ( 1 — “)  O')  — (7) — f7]  5 

rjo  =  c?(/>(1_“)[c2(,')_c"w_I/1  -  c0e“C([e“(/)“c3(/)“£/], 

rC0  =  (1  -  a)c?(/>(I-a)[c"W-c?«-£/l  +  occoe-^-^-u] 

(28) 


it  can  be  simplified  as  a  more  compact  matrix  form. 


Cl,/— i 


Cl, 


C1 ,7+1 


[1  0  0]  C2J- 1  +  [  —(2  +  h2<P2ro)  -h2$2rCo  h2&2rcJ  c2J  +[1  0  0]  C2J+1 


,c3j-l  , 


=  [~h202ro  - hz<PzrC0  hz0zr, 


(<j\ 

Cij 

VhJ 


,C3  J, 


,c3j+l  , 


+  h2<P2rj0.  0  <  j  <  N 


(29) 


In  the  same  way,  model  Eqs.  (18)  and  (19)  can  be  approximated  as  Eqs.  (30)  and  (31),  respectively. 


Cl,  7-1 


Cl, 


C1 ,7+1 


[0  1  0]  c2J-i  +[-h2(p-1rQ  —(2  +  h2<P2rC0)  h2&\rco]  c2J  +  [0  1  0]  c2j+i 


.  c3  / —  l 


,C3  J 


=  [  -h202ro  - h2<P2rC0  h2<P2rC0] 


f<j\ 
VhJ 


+  h2&\rj0,  0  <j  <N 


,c3j+l  , 


(30) 


Cl  ,7-1 


Cl, 


Cl,  7+1 


[0  0  1]  c2 j—  i  +  [h2<P2r0  h2<P2rCo  -(2  +  h2<P2rCo)  ]  c2j  +  [0  0  1]  c2j+  i 


.  c3  / —  l 


,C3  J 


,c3j+l  , 


=  [h2<t>\r0  h2<P2rC0  -h2$2rco} 


f<j\ 
VhJ 


-h2$22rj0,  0  <j<N, 


(31) 


28 


T.  Duan  et  al.  /  Journal  of  Power  Sources  107  (2002)  24-33 


Letting 


Eqs.  (33),  (35)  and  (36)  can  be  rewritten  as  a  matrix  form 


C(j)  =  (cij  C2j  c3J)  ,  0  <j<N  (32) 

then  Eqs.  (29)— (3 1 )  become 

A(j)C(j  -  1)  +B(j)C(j)  +  D(j)C(j  +  1)  =  G(j), 
0<j<N,  (33) 

in  which 


m 


m= 


m  = 


G(j)  = 


~h2$2rC0 
-(2  +  h2<I>\rC0) 


h2<t>\rc o 


h2$2rC0 
h2*\rC0 

-(2  +  /t2^rco) 


-hz<P~r0  -hz<P~rC0 
-h2$\ro  —h2<P2rCo 
h2022r() 

h2<P2- 


h2<P2rC(j  \  (c°i j\ 


h2<P2rC0 


h2*W 
~h202rco , 


-2j 

V%) 


+ 


rj0 
h2$\rj0 


(34) 


-h&in  <L 

Two-point  boundary  conditions  are  approximated  by 
three  point  forward  or  backward  difference  accurate  to 
0(h2)  and  are  also  expressed  as  a  compact  matrix  form 

fi(0)C(0)  +£>(0)C(1)  +XC{2)  =  G(0),  j  =  0  (35) 

and 

YC{N  -  2)  +  A(N)C{N  -  1)  +  B(N)C(N)  =  G(N ), 
j  =  N  (36) 

where 


m  = 


x  = 


Y  = 


'0  0  o' 

£>(0)  =  [  0  —4  0 

0  0  0. 


(37) 


3  0  0 
B(N)  =  |  0  3  0 
0  0  3 


(38) 


'B{ 0)  £>(0)  X 
A(l)  B(l)  D(  1) 

m  bo)  do) 

A(N  -  1)  B(N-l)  D(N  —  1) 

Y  A(N)  B(N) 


*  C( 0)  ■ 

■  G(0)  ■ 

C(l) 

G(  1) 

CO) 

= 

GQ) 

C{N-  1) 

G(N  —  1) 

.  C(N)  _ 

.  G(N)  _ 

The  block  matrixes  B(j)  and  vectors  G(/)  at  node  points 
0  <j<N  are  not  determined  due  to  the  variables  at  expanded 
points  (with  superscript  0)  being  unknown.  As  a  result,  the 
tridiagonal  matrix  algebraic  Eq.  (39)  cannot  be  solved  directly. 
But  if  these  unknown  values  at  expanded  points  are  guessed  as 
trial  values,  then  Eq.  (39)  can  be  solved.  For  a  tridiagonal 
matrix  algebraic  equation,  a  compact  forward  and  backward 
substitution  algorithm  can  solve  it  quickly.  This  algorithm,  at 
first,  decomposes  the  tridiagonal  matrix  into  a  LU  form.  Then 
in  the  forward  substituting,  an  intermediate  vector  will  be 
acquired,  and  in  the  backward  substituting,  the  variable  vector 
can  be  obtained.  In  all  the  process  of  LU  decomposition  and 
forward  and  backward  substitution,  the  Gauss-Jordan  full 
pivoting  algorithm  is  needed  to  inverse  the  block  matrixes. 
Then  the  block  matrix  in  the  lower  and  upper  matrix,  the 
intermediate  vector,  and,  at  last,  the  variable  vector  can  be 
solved.  These  solved  variables  are  based  on  guessed  values, 
they  are  not  the  real  solution  for  the  problem.  They  are  set  as 
new  guessed  values  in  the  block  matrixes  B(j)  and  G(j)  and 
reproduce  the  whole  procedure.  The  iteration  will  not  stop 
until  the  setting  accuracy  is  reached,  which  makes  the  solu¬ 
tions  satisfy  the  whole  difference  equations  at  each  node  point. 

The  procedure  is  illustrated  as  follows.  When  the  values 
for  determining  the  B(j)  and  G(/)  are  guessed,  then  the 
coefficient  matrix  in  Eq.  (39)  can  be  decomposed 

~B( 0)  £>(0)  X 
A(l)  fi(l)  D(l) 

A0)  B0)  DO) 

A(N  —  1)  B(N-  1)  D(N-  1) 

Y  A(N)  B(N) 
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G(  0) 
G(l) 


29 


m 

A(l)  5(1) 

m  m 

A(N  —  1)  B(N  —  1) 

Y(N)  A(N )  B(N)_ 

y 

L 

7  5(0)  X 
I  D(  1) 

/  £>(/)  (40) 


/  D(N-l) 


Through  the  shift  of  matrixes,  the  block  matrix  in  the  L  and 
U  matrixes  can  be  solved.  For  j  =  0,  1,  we  have 

5(0)  =  5(0),  X  =  5  ‘(0)X,  5>(0)  =  5  1  (0)D(0) 

(41) 

A(l)  =A(1),  5(1)=5(1)-A(1)D(1), 

/)(1)  =  5  '(1)5)(1)  —  5  ‘(1)A(1)Z  (42) 

As  the  same  way,  for  other  node  points  at/  =  2, ... , 
N-  1, 

A(j)  =  A(j),  B(j)  =  B(j)  -  A(j)D(j  -  1), 

D(j)  =  B~1(j)D(j)  (43) 

can  be  solved.  At  the  last  node  point  j  =  N,  we  obtain 

Y(N)  =  Y,  A(N)  =  A(N)  -  Y(N)D(N  -  2), 

B(N)  =  B(N)  -  A(N)D(N  -  1)  (44) 

Then  Eq.  (41)  can  be  solved  through  the  solution  of  two 
sets  of  matrix  algebraic  equations. 

-5(0) 

A(l)  5(1) 

m  m 

A(N  —  1)  5(51-1) 

Y(N)  A(N)  5(A)  _ 


"  W(0) 

W(  1) 

W(j)  =  G(j)  (45) 

W(N-  1)  G(N  —  1) 

W(N)  J  L  G(N) 

v  j  v  j 

W  G 

and 

"/  0(0)  X 
I  D(  1) 

/  m 


I  D(N  —  1) 


c  w 

The  intermediate  vector  W(j),  j  =  0. . . . ,  N  in  Eq.  (45)  can 
be  determined  by  forward  substituting 


be  determined  by  forward  substituting 

W(0)  =  5_1(0)G(0)  (47) 

W(l)  =  5  ‘(l)G(l)  -  5_1(1)A(1)1T(0)  (48) 

W(j)  =  B  l(j)G(j)  -  B  l{j)A(j)W(j  -  1)  (49) 

W(N)  =B~l(N)G(N )  —  B  1  (N)  Y  (N)  W(N  —  2) 

—  B^1  (N)A(N)W(N  —  1)  (50) 

With  backward  substituting,  the  solution  vector  C(f),  j  =  0, 
N  are  solved. 

C(N)  =  W{N )  (51) 

C(N  —  1)  =  W(N  —  1)  —  D(N  —  l)C(N)  (52) 

C(j)  =  W(j)  —  D(J)C(J  +  1)  (53) 

C(0)  =  W(0)  -  5>(0)C(1)  -  XC{ 2)  (54) 
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The  solved  C(j)  are  as  new  guessed  values  to  determine 
the  block  matrixes  B(j)  and  G(j ),  j  —  1 .....  A'  —  1  and  to 
continually  carry  out  the  iteration,  until  the  approximated 
solutions  are  obtained. 


4.  Results  and  discussion 

4.1.  Convergence 

To  test  the  solving  technique,  we  chose  a  set  of  physical- 
chemical  parameters  shown  in  Table  1  that  are  common  used 
for  hydrogen  oxidation  in  a  fuel  cell  anode,  then  the  follow¬ 
ing  dimensionless  model  parameters  are  obtained 

c0  =  0.5714,  a  =  0.5,  U  =  0, 

<Z>2  =  0.1406,  <Z>2  =  <|>2  =  0.8277 

Different  currents  density  (0.2-1. 6  A  cm-2)  at  two-point 
boundary  that  make  kt  and  k2  change  between  0.2191- 
1.7526  were  adopted.  A  Fortran  95  code  was  developed 
to  do  the  numerical  solution  on  a  micro  computer  (Dell, 
Dimension  4100,  Pentium  III,  260  MB  memory).  After  the 
number  of  node  points  was  determined  and  a  set  of  proper 
initial  values  were  guessed  for  each  variables  at  all  the  node 
points,  the  iterating  process  was  carried  out.  It  converged 
very  quickly  in  just  two  loops.  The  numerical  results  are 
illustrated  as  Figs.  2-5. 

It  was  found  that  a  large  number  of  node  points,  for 
example,  >1000,  is  required  to  converge  the  iteration.  We  set 
it  as  2000.  Setting  too  many  node  points  is  the  need  for  our 
assumption  in  the  Taylor  expansion  that  the  nonlinear 
difference  equation  at  node  point  j,  depends  only  on  vari¬ 
ables  at  the  three  points.  The  iteration  precision  was  chosen 
as  10~4.  The  zero  in  the  right-hand  side  of  the  sets  of 
nonlinear  algebraic  equations  in  the  form  of  Eq.  (26)  is 
approximated  by  this  value.  If  all  the  absolute  values  of 
computation  for  the  terms  in  the  left-hand  side  are  not 
>  1 0  4,  then  the  iteration  processes  stop.  For  the  case  of 
2000  node  points,  correspondingly,  there  are  6000  algebraic 


Table  1 

Common  physical-chemical  parameters  in  modeling  the  hydrogen 
oxidation  in  a  fuel  cell  anode 


Parameter 

Value 

a  (cm-1) 

250 

CH+  (mol  cm-  ) 

0.4  x  10~4 

C„2  (mol  cm-3) 

0.7  x  10~4;  0.38  x  10~4  for  <2>  =  0.15 

DHl  (cm2  s-1) 

2  x  10~4 

/Ceii  (A  cm-2) 

0.2-1. 6 

k°  (cm  s-1) 

0.45 

/  (cm) 

5  x  10~4 

T  (K) 

353.15 

U 

0 

a 

0.5 

k  (Q-1  cm-1) 

0.3  x  10_1 

o  (Q-1  cm-1) 

0.3  x  10_l 

5000  V 


Fig.  2.  Distribution  of  electrochemical  reaction  rate  of  hydrogen  along  x  at 
different  k\  (current). 

equations  like  Eq.  (26)  except  those  at  the  two-point  bound¬ 
ary.  Like  all  trial  and  error  iteration  algorithm,  different 
initial  estimations  yield  different  results,  in  order  to  make 
the  solution  meet  the  physical  meanings  in  the  system, 
following  current  criteria  are  introduced  in  our  computing 
to  determine  the  initial  trial  solution  and  check  the  results 


cT  0.85  4 


0.8  4 


Fig.  3.  Distribution  of  dimensionless  hydrogen  concentration  C\  along  x  at 
different  k\  (current). 
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3 


0  0 


Fig.  4.  Distribution  of  dimensionless  potential  in  matrix  along  x  at 
different  kx  (current). 


Fig.  5.  Distribution  of  dimensionless  potential  in  pore  solution  c3  along  x 
at  different  ky  (current). 


12  =  ~K 


d&2 

dX 


=  0 


x=o 


f'di2  , 

/cell  =  —6x 

,/n  dx 


(58) 

(59) 


For  Eqs.  (55)-(58),  the  five  point  formula  of  numerical 
differential  are  adopted  to  calculate  the  potential  gradient. 
For  Eq.  (59),  complex  Simposon’s  method  is  used  to  calcu¬ 
late  the  integration  along  the  catalyst  layer.  All  the  calcula¬ 
tions  in  these  five  equations  are  based  on  the  data  obtained 
from  the  numerical  solution  of  the  model  equations. 


4.2.  Characteristics  of  electrochemical 
reaction-diffusion 

The  characteristics  of  reaction-diffusion  problems  are 
usually  analyzed  by  the  Thiele  modulus  <P  L 10] .  We  changed 
the  model  parameters  that  is  shown  in  Table  2  to  get  different 
Thiele  modulus  <P  in  a  range  of  0.008-10,  and  corresponding 
electrochemical  Thiele  modulu  (I>t  and  <P2  in  a  range  of 
0.01-8.06  and  0.02-10.95,  respectively.  The  numerical 
solutions  converge  to  the  results  very  fast  too,  two  loops 
of  iteration  for  all  case.  The  results  are  showed  in  Figs.  6-9. 


Fig.  6.  Dependence  of  reaction  rate  on  Thiele  modulus,  the  parameters 
corresponding  to  these  values  of  <P  are  shown  in  Table  2. 


When  the  Thiele  modulus  is  low,  such  as  <P  =  0.008,  the 
resistance  of  diffusion  can  be  neglected,  the  distributions  of 
reaction  rate  and  hydrogen  concentration  are  nearly  uniform. 
When  the  Thiele  modulus  becomes  larger,  the  diffusion 


Table  2 

Parameters  for  studying  the  dependence  of  hydrogen  electro-oxidation  reaction-diffusion  on  Thiele  modulus  in  Figs.  6-9a 


0 

<t>i 

DHl  (cm2  s  ') 

k°  (cm  s  ') 

/  (cm) 

o  (Q  1  cm  L) 

k  (Q  1  cm  !) 

0.008 

0.012 

0.018 

2.0  x  10~4 

0.45 

0.1  x  10~4 

0.6  x  10~‘ 

0.3  x  10~‘ 

0.8 

0.41 

1.85 

0.1  x  10~4 

0.4  X  102 

0.25  x  10~4 

0.3  x  10~‘ 

0.15  x  10~2 

3 

1.41 

4.71 

0.08  x  10~4 

0.8  X  102 

0.6  x  10~4 

0.3  x  10~‘ 

0.268  x  10~2 

5 

1.73 

6.42 

0.08  x  10~4 

0.8  X  102 

1.0  x  10~4 

0.55  x  10~' 

0.40  x  10~2 

10 

8.06 

10.95 

0.08  x  10~4 

0.8  X  102 

2.0  x  10~4 

0.1  x  10~‘ 

0.55  x  10~2 

a  Non-changing  parameters: 

CH+  =  0.4  X  IQ" 

“4  mol  cm-3;  Ch2 

=  0.65  X  10~4  mol  cm" 

~3;  /cell  =  1.6  A  cm"2; 

T  =  353.15  K;  U  =  C 

);  a  =  0.5. 
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0  0.2  0.4  0.6  0.8  1 

Fig.  7.  Dependence  of  dimensionless  hydrogen  concentration  on  Thiele 
modulus,  the  parameters  corresponding  to  these  values  of  <P  are  shown  in 
Table  2. 


Fig.  8.  Dependence  of  dimensionless  potential  in  matrix  on  Thiele 
modulus,  the  parameters  corresponding  to  these  values  of  <P  are  shown  in 
Table  2. 


Fig.  9.  Dependence  of  dimensionless  potential  in  solution  on  Thiele 
modulus,  the  parameters  to  these  values  of  0  are  shown  in  Table  2. 


resistance  increases,  as  a  result,  the  reaction  is  not  uniform 
along  the  catalyst  layer,  the  larger  the  Thiele  modulus  is,  the 
more  non-uniform  the  reaction  rate  is.  When  the  Thiele 
modulus  is  equal  to  10,  the  reaction  rate  and  hydrogen 
concentration  decreases  rapidly  in  the  region  near  the  sur¬ 
face  of  porous  catalyst  layer.  It  should  be  pointed  out  that  in 
high  Thiele  modulus  region,  the  rate  of  mass  transfer  by 
diffusion  cannot  meet  the  rate  of  electrochemical  reaction, 
the  process  will  lose  its  stability  and  become  a  dissipative 
system.  As  the  same  case  in  the  stability  issues  in  chemical 
reaction  engineering,  our  numerical  computing  also  found 
sensitive  dependence  of  the  solution  results  on  some  para¬ 
meters  when  the  Thiele  modulus  is  high.  It  has  been  reported 
by  a  large  amount  of  literature  that  in  open  chemical  reaction 
system,  especially  coupled  nonlinear  chemical  reaction  and 
diffusion  occur  simultaneously,  many  complex  phenomena 
will  take  place,  such  as  multiple  steady  states,  unstable  states 
and  self-generated  sustained  oscillations.  Our  present  work 
just  makes  a  beginning  in  this  field,  further  attention  should 
be  paid  to  it  because  the  electrochemical  oxidation  reaction 
rate  of  hydrogen  is  very  fast  under  the  catalysis  of  platinum 
catalyst.  When  a  fuel  cell  stack  is  needed  to  scale  up,  it  is 
necessary  to  know  the  dissipative  behavior. 

5.  Conclusion 

A  numerical  solution  method  is  developed  for  coupled 
electrochemical  reaction-diffusion  equations.  Even  though 
the  convergence  of  this  technique  depends  mainly  on  the 
large  number  of  node  points,  computation  experiences  on  a 
microcomputer  with  double  precision  real  only  eight  word 
bites  shows  that  the  procedure  converges  very  fast.  After  the 
Taylor  series  is  expanded  and  the  corresponding  iterating 
procedure  is  constructed,  the  computation  process  is 
stable.  Arranging  the  elements  in  the  coefficient  matrix 
to  block  matrix  form,  the  whole  coefficient  matrix  is  easily 
decomposed  to  lower  and  upper  matrix,  and  the  compact 
forward  and  backward  substitution  algorithm  based  on  the 
shift  of  block  matrixes  with  Gauss-Jordan  full  pivoting 
method  can  perform  the  numerical  calculation  quickly. 
Local  convergence  depends  on  the  first  trial  solution, 
current  criteria  are  required  to  make  the  solution  converge 
to  the  correct  results.  It  is  suggested  by  the  model  solutions 
that  dissipative  behaviors  in  a  electrochemical  reaction- 
diffusion  system  might  occur  when  the  Thiele  modulus  is 
high. 
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